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Abstract 

The novel contribution of this paper relies in the proposal of a fully im- 
plicit numerical method designed for nonlinear degenerate parabolic equa- 
tions, in its convergence/stability analysis, and in the study of the related 
computational cost. In fact, due to the nonlinear nature of the underly- 
ing mathematical model, the use of a fixed point scheme is required and 
every step implies the solution of large, locally structured, linear systems. 
A special effort is devoted to the spectral analysis of the relevant matrices 
and to the design of appropriate iterative or multi-iterative solvers, with 
special attention to preconditioned Krylov methods and to multigrid pro- 
cedures: in particular we investigate the mutual benefit of combining in 
various ways suitable preconditioners with V-cyclc algorithms. Numerical 
experiments in one and two spatial dimensions for the validation of our 
multi-facet analysis complement this contribution. AMS SC: 65N12, 
65F10 (65N22, 15A18, 47B35) 



1 Introduction 

We consider a single equation of the form 

du , , . 

— = V.(D(u)Vu) 7 (1) 

where D{u) is a non- negative function. The the equation is parabolic and it 
called degenerate whenever D(u) vanishes for some values of u. For the conver- 
gence analysis of our numerical methods, we will require that D(u) is at least 
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diffcrcntiablc and that D'(u) is Lipschitz continuous, while the existe nce of so- 
lution s is guaranteed under the milder assumption of continuity (see Vazquez! 
l2007ft . 

The classical porous medium equation (where D(u) is restricted to be a 
power law) or its generalized form have important applications in many fields 
of science. Their name arise from the model of Darcy's flow of a gas through 



a porous medium (D(u) 



KU 



where 7 is the specific heat ratio), but classical 



applications range from ground- water flow (Boussinesq equation) to spatial pop- 
ulation dynamics (where crowding effects require nonlinear degenerate diffusion 
terms). Moreover models based on equation (fTJ) have been proposed as useful 
approximations of more complex models like thin films motion (when disregard- 
ing surface tension), water-oil mixtures in porous medium, boundary layer in 
fluid flow past an obstacle, magma in volcanoes, etc. More recently they have 
been studied as limits of kinetic particle models, with applications for the dif- 
fusion in semiconductors. Finally we mention that some contrast-enhancement 
filters for image processing like the one b y Perona-Malik are based on JT]). For 



Vazquez! ([2007I . Chap 2 and 21) and 



more details on the aplications, see e.g. 
the references therein. 

The present investigation is part of the search for suitable numerical tech- 
niques to integrate for long times nonlinear, possibly dege nerate, parabolic equa- 



tions a ppearing in models for monument degradation fsee lAregba Driollet et al 



( 2004 )) when chemical/micro-biological pollutants are taken into consideration. 
We wish to point out that the techniques developed here have applications that 
go beyond the aforementioned models. For example, again in the area of conser- 
vation of the cultural heritage, they could be ada pted to numerica lly investigate 
the more complete sulfation model describe d in jAli et all (|2007h and the con- 
solidation model presented in IClarelli et all <)2009h . So me applications in th e 
field of monument conservation have been presented in Semplice et all ( 2009f ) , 
where the mathematical tools developed in the present paper are employed for 
forecasting marble deterioration. Of course, in such a context, given the wide 
variety of artefacts, an important challenge is the combination of the approxi- 
mation scheme with the related linear algebra solvers, in presence of complicate 
geometries and griddings. 

In the literature, degenerate parabolic equations have been discretized mainly 
using explicit or semi-implicit methods, thus avoiding to solve the nonlinear 
equation arising from the elliptic operator. A remarkable class of methods arise 
directly from the so-called non-linear Chernoff formula iBrezis fe; Pazyl (|1972l ) 
for time advancement, coupling it w ith a spatial discret ization: for finite diffcr- 
en ces the latter st u dy wa s started in lBerger et and for finite elements 



in iMagenes et all (|1987T ). More recently, another class related to the relax- 
ation approximation emerged: such numerical procedures exploit high order 
non-oscillatory methods typical of the discretization of conservation laws and 
their convergence can be proved making use of s emigroup argument s similar to 



those relevant for proving the Chernoff formula ([Cavalli et alL\2001v . 
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In this paper we start from the Crandall-Liggett formula 



U(t n , x) - AtL D {U(t n ,x)) = U(t n -\x), (2) 

where time has been discretized with steps At = t n — t™ -1 , and —Ld{') denotes 
the elliptic operator u <— > —V • (D(ji)Vh). The computation of the numerical 
solution U (t n , x) requires to solve a nonlinear equation whose form is determined 
by the elliptic operator and the nonlinear function D(u), but the convergence i s 



guaranteed without restrictions on the time step At (jCrandall fc Liggettl Il97lf ) . 
Furthermore, due to the nonlinear nature of the underlying mathematical model, 
the use of a fixed point scheme is required and the choice of the Newton-like 
methods imp l ies th e solution at every step of large, locally structured (in the 



sense of iTiHil . 119981 ) linear systems. A special effort is devoted to the spectral 
analysis of the relevant ma trices and to the design of appropriate iterative or 
multi- iterative solvers (see Serra-Capizzanol 1993 ). with special attention to 



preconditioned Kry l ov methods and t o multigrid procedures ( see iGreenbaum 
(|l997l ) ; ISaadl (|2003h ; iHackbuschl (|l985f ) ; iTrottenberg et all (|200l[ ) and references 
therein for a general treatment of iterative solvers). Although most of the 
analysis is developed in the one-dimensional case (from Section [5] to Section 
U, we also indicate in Section [5] how to generalize our approach to two spatial 
dimensions, we also perform numerical experiments for the validation of our 
analysis in both settings (see Section @] and Section [SJ respectively). 

The paper is organized as follows. In Section [5] we couple the time dis- 
cretization @ with a spatial discretization based on finite differences and set 
up a Newton method for the resulting system of nonlinear equations. We re- 
port the explicit form of the Jacobian appearing in the Newton iterations and 
we prove the convergence of the Newton methods under a mild restriction on 
At. In Section [3] we consider various iterative methods for the solution of the 
inner linear systems involved in the Newton method. A brief spectral analy- 
sis of the related matrix structures is provided in order to give an appropriate 
motivation for the good behaviour of the proposed iterative and multi-iterative 
solvers. In SectionUwe perform some numerical tests. In Section[S]we describe 
a generalization of the previous methods to a two-dimensional case and perform 
numerical tests in this setting too. Finally, a conclusion section with a short 
plan for future investigations completes the paper. 



2 Numerical methods in one dimension 

In order to discretize equations like ([1]) , we will employ a time semi-discretization 
given by the Crandall-Liggett formula and a space discretization based on finite 
differences, explained in the following subsection. The latter numerical choice 
leads to a system of coupled nonlinear equations that need to be solved at each 
discrete timestep in order to compute the solution of the PDE: this is achieved 
using the Newton method, as detailed in Subsection 12.21 where we also prove 
and comment convergence results. 
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2.1 Finite difference discretization 



We take into consideration a standard discretization in space using finite dif- 
ferences. Denoting x$ = a + £h, we consider N + 2 points with equal spacing 
h = (b — a)/(N + 1) in the interval [a, b] and we denote by w£ the approximate 
solution at time t n and location Xk, where k = 0, . . . , iV+1. Let u" be the vector 
containing the collection of the unknown values u^. When no potential confusion 
can arise, we will sometimes drop in both notations the superscript indicating 
the time level. Of course, when considering Dirichlet boundary conditions, the 
values uo and un+i are known and can be eliminated by the equations, leaving 
a vector of unknowns u" that contains only uJJ for k = 1,...,N. Boundary 
conditions of Neumann or Robin type can be treated in similar ways. 

We choose a standard 3-points second order approximation of the differential 
operator (D(u)u x ) x . Denoting with the subscript £ the evaluation at the point 
X£, we have that: 



d_ 

dx 



D(u) j+1/2 m\j+l/2 _£) Wi-l/2 m\j-l/2 



du I 



D(u)j+i/2(u j+ i - Uj) - D(u) :) _ 1/2 (u j - Uj-x) 
h 2 



o(l) 



(D(u j+1 ) + D( Uj ))(u j+1 - Uj ) - (D(uj) + Djuj-^juj - Uj-i) 

2b? 



(3) 



where the o(l) error term is of order h 2 under the assumption that the compo- 
sition 

= £>(«(■)) 

is at least continuously differentiable, with Lipschitz continuous first derivative. 
Putting together all the contributions for different grid points, we end up with 
Z/£>( U )U, where the tridiagonal matrix 



— D\/2 - D 3 / 2 



D 



3/2 



conatins the values 



D3/2 

-D3/2 - D5/2 D 5 / 2 



D 



5/2 



Dn-1/2 
Dn-1/2 —D n _ 1 / 2 — D N+1 / 2 

(4) 



D 



i+1/2 



D{u j+1 ) + D{ Uj ) 



j = 0,...,N, 



and thus depends nonlinearly on the Uj's. It should be noticed that the latter 
is a second order approximation of (f>{xj + i/ 2 ) since m£ differs from u(t n ,Xk) 
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by 0(h 2 ) thanks to the second order scheme and since, by standard Taylor 
expansions, we have 

D{u j+ i) + D{uj) ^(^J + ^Xj) 2 
"3+1/2 = 2 = 2 ' ' 

h 2 

= <Kx j+1/2 ) + -j(l> xx {r){h,j)) + 0(h 2 ) 

= <f>(x j+ i/ 2 ) + 0(h 2 ), T)(h,j) € (xj_i,Xj), 

under the mild assumption that xa; (-) is a bounded function. Of course, the 
same conclusion holds if <f> x {') is Lipschitz continuous. 

In the following, we denote by tridiag fc [/3fc, ctfc, 7^] a square tridiagonal matrix 
of order N with entries /3 k on the lower diagonal, k = 2, • • • , iV, ot k on the main 
diagonal, k = 1, • • ■ , AT, and 7^ on the upper diagonal, fc = 1, • • • , N — 1. With 
this notation, L D(u) = tridiag fc [L> fc _i/ 2 , -D k _ 1/2 - D k+1/2 , D k+1/2 ]. We also 
denote with diag fe [a!fc] the square diagonal matrix with ctk on the k th row. 

As already observed L D ^ is a symmetric real tridiagonal matrix. Since -D(-) 
is a nonnegative function, the matrix —L D ^ is always positive semidefinite, 
beacuse it is weakly diagonally do minant by row, or equivale ntly thanks to the 
first Gerschgorin Theorem (see e.g. Golub fc Van Loan! 1 19961 ). Furthermore, we 



have positive definitcncss (i.e. invertibility) , at least for every N large enough, 
if in addition </>(•) has only isolated zeros in (a, b). In that case, for N large 
enough, the matrix is irreducible or block diagonal with irreducible blocks. In 
particular, when </>(•) is strictly positive in (a, b) then —Lo( u ) is positive definite 
and irreducible for any N . 

When introducing numerical methods for the approximation of the differen- 
tial equation we will encounter nonlinear systems involving the matrices — L D ( u y 
At that point more sophisticated (spectral) relations and features will be dis- 
cussed, when choosing the appropriate iterative solvers for the global linearised 
system (see Subsection 13. 1[) . For the moment we just observe that, thanks to 
the previous preliminary spectral analysis, all the classical iterative solvers like 
Jacobi and Gauss-Seidel (and their damped version with damping parameter 
belonging to (0, 2)) are all convergent for the solution of a linear system with 
such a coefficient matrix. The problem is that the spectral radii are very close to 
1, with a gap ranging between 0(N~ 2 ), reached by all these classical iterations 
with the only exception of the optimally damped Gauss-Seidel, and OjN -1 ) , 



reached for Gauss-Seidel with optimal damping parameter (see IVargal |1962|). 
When considering the whole system things become slightly better since the gap 
between the spectral radius and 1 reduces for all the considered procedures to 
0(N^ 1 ). However, as a partial conclusion, we can safely claim the considered 
iterations would be unacceptably slow and the search for specialised iterative 
solvers becomes mandatory. This latter is the main subject of Section [3l 
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2.2 The nonlinear system and the Newton iteration 

Following the Crandall-Liggett formula ©, in order to compute u" from u" _1 , 
we need to solve the nonlinear vector equation 



u 



At 



h 2 ^C(u") u 

and thus we set up Newton iterations for the vector function 
F(u) = u ^-L D{u) u - u"" 1 . 



(5) 



In the following, we denote u n ' s the s th Newton iterate for the computation of 
u". The generic partial derivative of F(u) is 



dF k At 
^7 <*i* -^ L 



\j,k 



At 
2h? 



$k-i,jD' k _i(uk-i - u k )+ 
+ h,j D 'k( u k-i - 2wfe + Ufc+i) + 
+ Sk+i,jD' k+1 (uk+i - u k ). 



so that the Jacobian is 



F'(u) =X N {u)+Y N {u), 
At 

-Ajv(u) = In — 7jiD( u )i 



Y N (u)=-^T N (u)d\ag k [D' k ], 



T N (u) = tridiag fe [u fe _i -u k ,u k ^ 1 - 2u k + u k+ll u k+1 - u k }- 



(6) 

(7) 
(8) 

(9) 
(10) 



The matrix X^(u) is symmetric positive definite and X m i n (X jy (u)) > 1, where 
Amin(^l) denotes the minimum eigenvalue of the matrix A. We note that the 
inequality is strict under assumption of isolated zeros. If diag fc [D' k ] is positive 
semidefinite, i.e. if D(-) is a smooth nondecreasing function, then, setting E 2 = 
diag fc E is a positive semidefinite diagonal matrix and Yjv(u) is similar to 

At Z7 

' 2h 2 ' 



-MeT n (u)E. Moreover, defining 



Yn(u) = -^f N (u)d\ag k [D k ], 



(11) 



Tjv(u) = tridiag fc [ti fc _i - u kl 0,u k+1 - u k ] = T N (u) - diag fc [u fc _i - 2u k + u k +i], 

(12) 

we have that Yat(u) is similar to ~^^ET^(u)E, with the latter being anti- 
symmetric, which implies a pure imaginary spectrum. 

In the following we will denote by || • || the Euclidean norm for vectors and 
the induced spectral norm for matrices. 

Remark 2.1. If u is a sampling of a solution u of ^ at least continuous and 
uj u (-) denotes its modulus of continuity, then 



\Y N (u)\\ < e u {At,h) 
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with 

e u {At,h) = A-^\\D'{u)\\ 00 uj u {h). 

In order to deduce the latter, it is enough to recall that for normal matrices 
the spectral norm is bounded by any induced norm and in particular by the the 
matrix norm induced by the infinity vector norm. In particular, if u is Holder 
continuous with exponent a e (0, 1] and constant M > 0, then the estimate 
above can be written as 

e u {At,h)=AM-^-\\D>{u)\\ 00 . 

Of course ifu is continuously differentiable, then we find a = 1 and M = ||u'||oo- 
Furthermore, when u is two times continuously differentiable, we notice that 
Uk-i — 2uk + Mfe+i = h 2 u"(£k) which leads to a more refined expression i.e. 



e u (At, h) = ^^Hu'lloo + At||u"||oo) \\D' 



Wile 



Finally, if we are interested in evaluating ||Yjv(u)|| where u is an approximation 
to the true solution u (this happens naturally in the numerical process discussed 
in the present section), then 

||Yat(u)|| < eu(At, h) < e u {At, /i)+4^||D / (u)|| 00 ||u-«|| o = 4^||£>'(u)|| (X) (u v (h) + \\u- u||oc) . 

Hence, since we are using second order formulae, the error \\u — = 0(h 2 ) 
and therefore ||Yjv(u)|| is dominated byuj u (h), which is of order h if the solution 
is Lipschitz continuous, that is 

\\Y N (u)\\ < 4M^||D / (u)|| 00 + O(At). 

In conclusion, we can safely claim that the global spectrum of the Jacobian F'(u) 
is decided, up to small perturbations, by the matrix Xn(u). For making more 
explicit the latter statement, if we assume that At — Ch, where C > is 
independent ofh, then X m i n (X n (u)) > 1, ||Xjv(u)|| = 0(/i _1 ), while ||Yjv(u)|| = 
0(1). 

In order to prove the convergence of the Newton method, we first consider 
some auxiliary results. 

Lemma 2.1. For a generic matrix A, the minimum singular value is 

( A + A T \ 

tfmin(A) > A min . (13) 



Proof. Consider the symmetric matrix 

A 1 ' 



B = 



A 
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with eigenvalues Xi{B) = <J\{A) > ••• > \n{B) = <tn{A) > Xn+i(B) = 
—<7n{A) > •■• > \2n{B) = —<t i (A), since the Schur decomposition of B is 
easily written in terms of the singular value decomposition of A (see iBhatial 
19971 iGolub fc Van Loanl . [l996h . Let V be a vector space and x e R N such that 



| >c| [ 2 > 0. Thanks to the minimax principle (jBhatial . 119971 ). we obtain 



/ a\ . , m . y T By . y T By . x T (A T + A)x f A + A 1 

AA) = An(B) = max mm — = — > mm — ^ — = mm — -= = A m i n — 

dim(V)=W yeV y 1 y y=[x. T , x T ] T Y y xeK™ 2x J x \ 2 



xt 



□ 

Remark 2.2. The proof technique used for bounding from below the minimal 
singular value of a matrix A is part of a more general framework useful for 
refining, when necessary, the estimates. In fact, in general it can be proved that 
for any complex-valued matrix A the minimal singular value is not less than the 
distance d r of any straight line r separating the numerical range of A from the 
complex zero. Therefore a better estimate can be obtained by computing the sup 
(that we call d) of d r , over all straight lines that induce the separation. In our 
case we used the fact that the real part of A (that is Ke(A) — (A + A T )/2) is 
positive definite and so our straight line becomes the set of all complex numbers 
having real part equal to A m i n (Re(A)). The estimate could be poor since the latter 
straight line is not necessarily tang ent to the num erical range (a convex set by 
the Toeplitz-Hausdorff theorem, see Bhatici 1 199$ ) ): thus d could be much larger 



than d r . However in our setting such an estimate is already very satisfactory, 
as also stressed by the numerical experiments. 

Proposition 2.3. Consider F(u) as defined in ([5|) ; where u is a sampling 
(at a given time t) of a solution u of ^ with D differentiable and having first 
derivative Lipschitz continuous. If, in addition, u is differentiable with Lipschitz 
continuous first derivative, then 

H^'(u)- 1 1| < l + 0(Ai). (14) 

When using the induced l°° norm, we have 

llFvri^d (15) 

for h sufficiently small and under the additional assumption that At < Cooh for 
some Coo > 0. 

Proof. For the sake of notational simplicity, we set A = F'(u). First of all we 
write the symmetric part of A as 

A + A T 

= X N (u) + Z N (u) , 

where 

Z N {u) = -— ^ tridiag fe [{D' k _ x - D' k ){u k -i - u k ), 2D' k {u k -x - 2u k + u k+1 ), {D' k+1 - D' k ){u k+1 - u k )] 
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By the regularity of D and u, we have that every entry of Zn(u) is of order jph? 
that is 0(At) and hence ||^jv(u)|| = O(At). Thus, recalling that A min (Xjv) > 1, 
it holds 

'A + A 1 "- 



> 1 - CAt 



(16) 



for some C > 0, that contains the infinity norms of the first derivatives of u and 
D and their Lipschitz constants. Using Lemma I2TT1 



Inequality (TT^|) now follows combining (jTTJ) and (fTf 

For the proof of the the estimate in l°° norm, we note that 



(17) 



F'(u) =X N (u)- ^tridiag fc 



u'(i k ),0(h),u'(i k ) diag fc 



where € [^.iit], e [x^, a^+i] and the constant in the 0(h) contains the 
Lipschitz constant of u'. We split F'(u) as 



F'(u) = -l(Z N - W n ) , 



where 

Z N = diag fc [z fe ] = diag fc 
W N = tridiag fc 
From ([18]), we have 



^+D k _ 1/2 + D k+1/2 +0(h 2 )D' k 
Dk-1/2 + \u'£ k )D' k _ x ,Q,D k+l/2 + \u'{l k )D' k+1 



(18) 



[f'^t 1 = -(i - z-^WnY 1 Z N \ 



For the factor Z N , it holds 



V AT II oo 



+ D k+1/2 + f t (1 + 0(At)£>£) 



< c 



At 



(19) 



(20) 



for /i sufficiently small and assuming that At < Cooh, recalling that D(-) > 0. 
For the remaining factor (/ — Z^Wn) -1 , we note that 



Z N Y W N = tridiag fc 



D k - 1/2 + H D 'k + Q(h 2 ) D k+ i/2 + H D 'k + Q(h 2 ) 

1 



Zk 



Zk 



and hence 

H-Zjv^Wivlloo < max 



D k -i /2 + D k+1/2 + 0(h) 



D k - 1/2 + D k+1/2 + ^-h(l + 0(At)D' k ) 



< a < 1, 
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for Coo > sufficiently small. Thus the spectral radius of Z 1 W is p(Z 1 W) < 
1 and we have 

{I - ZJ+WnY 1 = Y,{Z]j L W N ) S => \\ (I - Z^WnY 1 Wo. < — -. 

3=0 

(21) 

Finally, combining (gl) and (gD} with ([X5|). the {T5J holds with C\ = j^. □ 

Remark 2.4. TTie above result, with minor changes, can be proved under weaker 
assumptions. Indeed if both u(-) and D(-) are continuously differentiable, then 
every entry of Zn(u) is of order 



O ^max j— u> U '{h), -^Wd'(/i)II u 'I|oo j J 



wii/i denoting the modulus of continuity of a given function v. Therefore 

with the choice of At proportional to h and setting a(h) = max{uv(/i), u>D>{h)} = 
o(l), we find 

1 |T 1 >1-Ca(h) 

and by Lemma \2.1\ ||v4 _1 || < 1 + Ca(h). Furthermore, if we require that u is 
only Lipschitz continuous then the inequality regarding the norm of A^ 1 reads 
as ||^4 _1 || < C where C linearly depends on the Lipschitz constant of u. Finally, 
the same results can be obtained with minor changes, when using the induced 
l°° norm. 

Remark 2.5. In general, the solution u to ^ is not smooth, but only piece- 
wise smooth with a finite number of cusps. For instance with D{u) = u rn and 
continuous data with piecewise continuous derivative, the derivative of u is not 
defined in a finit e number of po ints in ID and in a finite number of smooth 
curves in 2D; see I Vdzauea 12001). The latter implies that the related matrices 
have the same features up to low rank correction terms whose cumulative rank 
is 0(N d ~ 1 ) if the equation is in d dimensions. 

Since the Crandall-LigKett formula d oes not induce any restriction on the the 



timcstep At (jCrandall fc Liggetu 119711 1 , we have only to prove the convergence 



of the Newton method. We are interested in the choice At = Ch for a constant 
C independent of h, which gives a method which is overall first order convergent. 
This is no restriction due to the presence of singularities at degenerate points: 
higher order methods would be more computationally intensive without reaching 
their convergence rate, even if in practice a certain reduction of the error is 
expected. 

Indeed, concerning the stopping criterion |ju" +1 ' s+1 — u™ +1 ' s || < e for the 
Newton method, the following observation is of interest. Since the method is of 
first order in time At can be chosen equal to h, it is sufficient to set e = c ■ h 
where c is moderately small constant independent of h. In fact, more precision 
will be useless in practice and would make the Newton process more expensive 
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by increasing the iteration count. The following result is a classical tool (see 
Ortega fc Rheinboldt dl970h ) for handling the global convergence of the Newton 
procedure. 

Theorem 2.6 (Kantorovich). Consider the Newton method for approximating 
the zero of a vector function F(u), starting from the initial approximation U^°K 
Under the assumptions that 



F'(u< 0) ) 
>( U (°)) 



< 



^(u (0) )|| < v , 



and that 



|F'(u)-F'(v)||< 7 ||u-v|| 



Pvi < 2 . 



(22a) 

(22b) 
(22c) 

(23) 



the method is convergent and, in addition, the stationary point of the iterations 
lies in the ball with centre and radius 



i - yi - 2/^7 

Pi 

For the choice At — Ch we can prove the following result. 

Theorem 2.7. The Newton method for F(u) defined in ([5]) for computing u" 
is convergent when initialised with the solution at the previous timestep (i.e. 
u n, ° = u n ~ 1 ) and for At < Ch, for a positive constant C independent of h. 



Proof. We will make use of the Kantorovich Theorem 12.61 so we need the esti- 
mates l|22p and to show that (f2"3")l is satisfied. We will use the V vector norm 
ll u llp = 12 \ v j\ P an d the related induced matrix norms. When p = 2 we find 
the Euclidean vector norm and the induced spectral norm; in general they are 
simply denoted as || ■ ||. 

Concerning (|22a[) . Proposition 12.31 and the assumption At < Cooh imply 



/3<Ci, 

at least for p — 2, oo, C\ — C\(p) and h small enough. 
Regarding (|22b|) 



(24) 



[^'(u"- 1 )] 1 J F(u"- 1 ) </3\\F(u 



At 



1 )' 



P u 



n-2 



for a constant C2 = C2(p) independent of h. The first equality in the previous 
calculation follows from ([5]) , while the second one is a consequence of the fact 
that u ra_1 is the stationary point of the Newton iteration for the previous time 
step and thus it satisfies 



u 



D(u™ 



-nil = U 
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It follows that 

i] = C 2 f3Ath- 1/p . (25) 

From now on wc consider only the || • ||oo norm i.e. p — oo, which leads to 
the most convenient estimate in ([23)1 and hence to the weakest constraint on 
the timestep At. 

For the Lipschitz constant of F', i.e., for estimating (|22c[) . observe that 
F'(u) — F'(y) is a tridiagonal matrix with two contributions: 



At 



F'(u) - F'(v) = -^(L D(n) - L D{V) ) + (Y N (u) - Fjv(v)), 
with Y/v(-) as in ([§]). The first term can be estimated as follows: 
||Ad(u) - Ad(v)||oo < 4||D'|| 0O ||u- v||oo. 



(26) 



(27) 



In order to check that the last inequality is satisfied, one observes that the sum 
of the absolute values of the entries in each row of L D ^ — Ljjm is smaller than 
the sum of 4 terms of the form 



|A=±i/2(u) - D k±1/2 (w)\ = 



D 



f Uk±i + u k 
V 2 



- D 



( Vk±l + Vk 

V 2 



\D'(C)\ < ||D'U|u- v|| 



For the second term in (f2l)l) , we have 



At 



\Y N (u)-Y N (v)\\ O0 <-^\\D'\ 



\M\ 



(28) 



where 



M = tridiag fc 



and hence 



(uk-i - u k ) - (vk-i - v k ) 

(uk-i - 2u k + Ufc+i) - (vk-i - 2w fc + Vk+i) 

{uk+i - u k ) - (vk+i - v k ) 

Halloo < 8||U- vHoo. 



(29) 

Replacing equation (gH]) in and combining (J3HJ) and ^7J) with (|2l)]). we 
obtain 

7<8||£>'|U^. (30) 
5]), and (1501 , Theorem 12.61 implies that 



Finally, combining equations 
Newton converges provided that 

l>C%C 2 8\\&U^>f}rn, 

i.e., At < Ch, for h sufficiently small and for C = min{C' 0O , l/(4Ci x /C 2 ||£>'||oo||£''||oo )} 
(essentially) independent on h. □ 
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3 Algorithms for the resulting linear systems 



At each Newton iteration, we need to solve a linear system whose coefficient 
matrix is represented by the Jacobian F'(u) with entries as in ([6]). In principle, 
the Jacobian is recomputed at each Newton iteration, so we are interested in 
efficient iterative methods for solving the related linear system. 

We note in passing that the form of the Jacobian matrix used here is very 
similar to the one that is obtained discretizing in space with P\ conforming 
finite elements. Thus the methods considered here can be to some extent gen- 
eralized to finite elements approximations. In particular, when considering real 
2D and 3D cases, the structure of the relevant matrices will depend heavily on 
the geometry of the domain, on the triangulation/gridding (often generated au- 
tomatically), and on the type of finite elements (higher order or non Lagrangian 
etc.). Therefore fast methods that are based on a rigid algebraic structure (e.g. 
of Toeplitz type) cannot be adapted because the structure is lost, in the gen- 
eral framework. However there exists a kind of information depending only on 
the continuous operator and which is inherited virtually unchanged in both fi- 
nite differences and finite elements, provided that the grids are quasi-uniform in 
finite differences and the angles are not degenerating in finite elements. Such in - 
formation consists in the locally Toeplitz structure (see ISerra-Capizzano (2006); 
Tilli (1998)) and in the related spectral features (conditioning, subspaces related 



to small eigenvalues etc.). We remind that these spectral features are conve- 
niently used when defining ad hoc preconditioned Krylov methods or multigrid 
algorithms, working uniformly well in one or more dimensions. 

In order to choose appropriate iterative methods for solving the jacobian 
linear system, we first analyse the spectral properties of the matrix F' (u) . This 
will lead us to consider preconditioned Krylov methods, multigrid and their 
combinations. 



3.1 Spectral analysis for the resulting matrix-sequences 

We start by introducing the notion of spectral distribution for a matrix sequence. 
Then we will briefly report a concise analysis of some delicate spectral features 
of the matrices involved in the definition of the Jacobian. Since the emphasis 
of this work relies in the computational aspects, we will not report all possible 
details, nuances, and generalisations of the spectral analysis. 

Definition 3.1. Let Co(R^) be the set of continuous functions with bounded 
support defined over the nonnegative real numbers, d a positive integer, and 9 
a complex-valued measurable function defined on a set G C M. d of finite ana 
positive Lebesgue measure fJ,(G). Here G will be often equal to (—ir,n) d so that 
e iG _ Y d with i 2 = — 1 and T denoting the complex unit circle. A matrix 
sequence {A/v} is said to be distributed (in the sense of the eigenvalues) as the 
pair (9, G), or to have the eigenvalue distribution function 9 ({A^} ~> (9, G) ), 
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if, VF G Co(C), i/ie following limit relation holds 

i N i r 

A m «JVp (Aj( ^ )) = MG)i fW))dt ' t = (ti,..., (31) 

Furthermore, a matrix sequence {An} is said to be distributed (in the sense of 
the singular values) as the pair (9,G), or to have the distribution function 9 
({An} ^ a (9,G)), if, VF e Co(Rq), the following limit relation holds 

1 N 1 r 

^onJ2 F ^ An ^ = -^G) //d^)!)^ * = (*i,---,*d)- (32) 

Along with the distribution in the sense of singular values/eigenvalues (weak*- 
convergence) , for the practical convergence analysis of iterative solvers we are 
also interested in a further asymptotic property called here the clustering. 

Definition 3.2. A matrix sequence {An} is strongly clustered at s e C (in the 

eigenvalue sense), if for any e > the number of the eigenvalues of An off the 
disk 

D(s,s) := {z:\z- s\< e} 

can be bounded by a pure constant q £ possibly depending on e, but not on n. In 
other words 

q £ (n,s) := #{\j(A N ) : Xj <£ D(s,e)} = 0(1), n -> oo. 

If every An has only real eigenvalues (at least for all n large enough), then s 
is real and the disk D(s,e) reduces to the interval (s — e,s + e). Furthermore, 
{An} is strongly clustered at a nonempty closed set S C C (in the eigenvalue 
sense) if for any e > 

q e (n,S):=#{\ j (A N ):\ j £D(S,e):=U se sD(s,e)}=0(l), n -> oo, (33) 

D(S,e) is the e -neighbourhood of S, and if every A N has only real eigenvalues, 
then S has to be a nonempty closed subset o/R. Finally, the term "strongly" is 
replaced by "weakly", if 

q e (n,s) = o(n), (q s (n, S) = o(n)) , n -> oo, 

in the case of a point s (a closed set S ), respectively. The extension of the notion 
in the singular value sense is trivial and is not reported in detail. 

Remark 3.3. It is clear that {A N } ~a {0,G) ({A n } ~ ct (9,G)) with 9 = s a 
constant function is equivalent to {An} being weakly clustered in the eigenvalues 
sense at s e C (in the singular value sense at s e Eq ). 

Now we briefly use the above concepts in our specific setting. Given the 
linear restriction on At imposed by the convergence of the Newton method 
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fTheorem l2.7[) , we are interested in the choice At = Ch for C > independent 
of h. However, for notational simplicity, here we assume At = h and note that 
analogous results hold for C > 0. 

Taking into account At = h and the re-scaling An = hF'(u), we consider 
the sequence {An} such that 



An — — i_D( u ) + Rn{u) 
R N (u) = hI N -~T N (u)d\ag k (D k ) 



with Tjv defined as in ([TT 

We have the following results, which are of crucial interest in the choice, in 
the design, and in the analysis of efficient solvers for the involved linear systems. 

Remark 3.4. The conditioning in spectral norm of An is of order N : this 
is implied dir ectly by Propos ition \2.3[ More in detail, by using the Bendixson 
Theorem ( see Stoer & Bulirscf\ , \2002\ , Theorem 3.6.1) the eigenvalues of An 



are localised in a rectangle having real part in [ch, C] and imaginary part in 
[— dh, dh] for some positive constants c, d, C independent of N. This statement 
is again implied by the analysis provided in Proposition \2.3\ for the real part, 
while for the imaginary part we note that (A — A T )/2 — — jtr\d\ag k [(D' k _ 1 + 
D£)(u fc _i-u fc ), 0, (D' k+1 +D' k )(u k+1 -u k )]. 

Remark 3.5. {A N } ~a, ct (0,G) with 9(x,s) = D(u(x))(2 - 2cos(s)), G = 
[a,b] x [0,2-7r] (distribution o f the zero o rder main term). The distribution of 
{Ld(u)} * s already known fsee \Till1 , \l99ti ), if we assume thatu is a sampling of a 
given function over a uniform grid. In our case the entries ofu represent an ap- 
proximation in infinity norm of the true solution, the latter being implied by the 
convergence of the method, and therefore by standard perturbation arguments we 
deduce {L D ( u j} ^a.ct (-9, G) with 9 and G as above. Moreover the tra c e nor m 
(sum of all singular values i.e. Schatten p norm with p — 1; see \Bhatid \l99i) ) 



of the remaining part Rn{u) is bounded by a pure constant C independent of N , 
when assuming that D' is bounded and u is at least Lipschitz continuous. The 
latter implies that the distribution o f {An\ is decided only by th e symm etric 
part that is, essentially, {Ld^} lsee \Golinskii & Serra-Cavizzanc . 200't , The- 



orem 3.4). Moreover any real interval containin g the spectrum \Lr>(„\} is also 



a str ong eigenvalue clustering set for {An} fsee \Golinskii & Serra-Capizzanc 
200± Corollary 3.3 and Theorem 3.5). 



Concerning the negligible term, we have that {Rn(u)} ^a,<t (OiG) and 
{R N (u)/h} ~ A)(T 0,G) with^{x,s) = l-D'(u(x))2isin(s)), G = [a,b] x [0, 2tt] 
(distribution of the first order term). 



Remark 3.6. Setting 



Pn — -iflfu) + hlN, 



we have {P^An} ~a,ct (equivalent, as already observed, to a weak eigen- 

value/singular value clustering): it follows from t he property of algebra o f the 
Generalized Locally Toeplitz (GLT) sequences (see \ Serra- Capizzana . 200(1 ). 
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In fact the preconditioned sequence {P^An} is also strongly clustered at 1 
both in the eigenvalue and singular value sense: we remark that the strong clus- 
tering property can be recovered via l ocal domain analy s is, by employing the same 
tools and the same procedure as in \Bertaccini et al. \200d. Theorem 3.7); see 



also S ection 3.1 and the conclusion section in Beckermann & Serra-Cavizzano] 



\'200j) and references therein. 

More in detail, by the Bendixson Theorem the eigenvalues of P^ 1 An are 
localised in a rectangle having real part in [1 — C\h, 1 + c^h] and imaginary 
part in [— d, d] for some positive constants C\, C2, d independent of N . This 
statement follows by noting that the eigenvalues of P^An belong to the field of 
value of Pjq 1 ^ 2 AjyPjy 1 ^ 2 . Considering a = x H P^ 1 ^ 2 A^P N 1 ^ 2 'x., for all x 6 C n , 
|x|| = 1, it holds that the real part of a is x H P n 1 ^ 2 (An + A N )P N 1 ^ 2 x/2 which 
belongs to [1 — cih, l + c 2 /i] by the analysis provided in Proposition \2.3\ A similar 
analysis stands for the imaginary part of a similarly to Remark\3.4\ 



Remark l3.6l is very important in practice, since it is crucial for deducing that 
the number of iterations of preconditioned GMRES is bounded by a constant 
depending on the precision, but not on the mesh that is on h (optimality of the 
method). This will be discussed in the next section. 

3.2 Iterative methods for the linear system 

In this section we consider some iterative methods for solving the linear system 
at each Newton step and study their convergence properties on t he matrix se - 
quence {^4 at}. A classical reference for the results quoted below is ISaadl (|2003h . 



GMRES We first consider the GMRES algorithm, since the antisymmetric 
part of An is negligible but not zero. 

Assume that An is diagonalisable and let An — \YAW~ 1 , where A = 
diag fe (Afe) is the diagonal matrix of the eigenvalues. Define 

e(™ 1 ) = min max |p(Afc)|. 
peP m :p(o)=i k=l,...,N 

Denoting with r^ the residual at the m th step of GMRES, it is a classical 
result that 

||r (m) || 2 < K2 (iy)e( m )||r(°)|| 2 . 

Thanks to Remark [3Tl k 2 (T40 ~ 1< Thus the GMRES convergence is deter- 
mined by the factor e^ m \ Thanks to Remark l3.4i it is possible to construct an 
ellipse properly containing the spectrum of An and avoiding the complex 0, so 
that when one applies GMRES to the matrix An, it holds that 

e (m) < (l-CV/T) m (34) 

for a positive constant C that is independent of the problem size N. 
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Similarly, using Pn as preconditioner, Remark 13.61 implies that 

e {m) < C m (35) 

for some C € (0, 1), independent of the problem size N. Even if the solution 
u is not enough regular to assure that the spectrum of PZ An belongs to 
[1 — c\h, 1 + c 2 h] x \\—d, d], the strong cluster at 1 leads in practice to the 
super-linear convergence. 

Conjugate gradient (CG) Let Sn = (An + A N )/2 be the symmetric part 
of A N and define ||x|| Sjv = || (5 JV ) 1/2 x|| 2 . Denoting k 2 (S n ) = \\S N \\ 2 \\ S N %, we 
recall the following classical result about the convergence of the CG: 



l|x m -x,|| Sw <2( ^g4 I ) W*o-x4s N , (36) 
\V^{S N ) + 1/ 

where x m is the approximate solution obtained at the m th step of the CG 
algorithm and x* the exact solution. 

Thus, combining (|3l)]) with Remark 13.41 we expect the CG algorithm to 
converge in 0(-/~N) iterations when applied to Sn- On the other hand, using 
Pjv as preconditioner, (|36[) together with Remark 13.61 imply that CG converges 
in a constant number of iterations, independently on the size N of the problem. 

Finally, according to Remarks 13.41 and 13.51 the antisymmetric part of An 
is negligible. Thus in practice one may apply the CG algorithm to the matrix 
An, expecting a convergence behaviour similar to that for Sn, in both the 
unpreconditioncd and preconditioned cases. 

Multigrid method (MGM) From Remark 13.51 we have that An has the 
same spectral behaviour of — L D ( u y Hence, if an iterative method is effective 
for L fl ( u j and robust, it should be eff ective also for An . Thi s is the case of 



MGM largely used with elliptic PDEs (jTrottenberg et all 120011 ) 



MGM has essentially two degrees of indetcrmination: the choice of the grid 
transfer operators and the choice of the smoother (pre- and post-smoother, if 
necessary). In particular, let -Pj l +1 be the prolongation operator from a coarse 
grid t + 1 to a finer grid i. We consider a Galerkin strategy: the restriction 
operator is (P£ +1 ) T and the coefficient matrix of the coarse problem is Ai + i — 
(Pl +l ) T AiPl +l , where Ai is the coefficient matrix on the i th grid. 

For the prolongation we consider the classical linear interpolation. We note 
that it is not necessary to resort to more sophisticated grid transfer operators 
since An is spectrally distributed as — L D r u y The restriction is the full-weight 
since, according to the Galerkin approach, it is the transpose of the linear inter- 
polation. Concerning the smoother damped Jacobi, damped Gauss-Seidel and 
red-black Gauss-Seidel are considered. 

Remark 3.7. The robustness of MGM could be improved in several way. A 
possibility is to use as post-smoother a damped method that reduces the error 
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in the middle frequencies whose could be not well dealt with the pre-smoother 
and the coarse grid correcti on. This is called as "intermediate iteration" in 
the multi-iterative methods i Serra-Cavizzanol , \l993l ). Another degree of free- 



dom is the number of smoothing i t eratio ns depending on the grid i. Indeed in 



Serra-Cavizzano & Tablino-Possi<\ {200 A) it is shown that a polynomial growth 



with i does not affect the global cost, that remains linear for banded structures, 
only changing the constants involved in the big O. 

In our present setting is not necessary to resort to the strategies described in 
the previous remark. In fact the method that achieves the smallest theoretical 
cost and that minimises the CPU times, for reaching the solution with a preas- 
signed accuracy e, is the simplest V-cycle with only one step of damped Jacobi 
as pre-smoother. The reason of the observed behaviour relies in the spectral fea- 
tures of our linear algebra problem: indeed, An can be viewed, after re-scaling, 
as a regularised weighted Laplacian since in the coefficient matrix one adds h 
times the identity (see the previous subsection). In this way the conditioning is 
not growing as N 2 as in the standard Laplacian but grows only linearly with N 
(see Remark GO}. 

Therefore the basic V-cycle, with one single step of damped Jacobi as pre- 
smoother, is alread y optimal for An, i.e. the number of iterations is independent 



of the system size (jTrottenberg et all 120011) . Moreover, as we will see in the 



numerical tests of section 14.21 the number of iterations for reaching a given 
accuracy is already very moderate. Therefore the additional cost per iteration, 
that should be paid for increasing the number of smoothing steps and for the 
use of a post-smoother, can not be compensated by a remarkable reduction of 
the iteration count. 

Finally, we stress that a robust and effective strategy is to use a multigrid 
iteration as preconditioner for GMRES as confirmed in the numerical experi- 
ments. In fact we showed that Pn is an optimal preconditioner for GMRES and 
the MGM is an optimal solver for a linear system with matrix Pn . 



4 Numerical tests 

In this section we consider as a test case the porous medium equation written 
in the form 

with homogeneous Dirichlet boundary conditions. Here m > 1, with m = 1 
corresponding to the heat equation. In particular we consider the exact self- 
similar solution 



u(t,x) = t~ a 



1-0- 



2 m 



{\x\t-°) 



m+l 



(38) 



due to Barenblatt and Pattle Vazquez! ( 2007 ). (The subscript + denotes the 
positive part). The experiments are carried out in Matlab 7.5. 
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Figure 1: h error at final time for N = 32, 64, . . . , 2048, m = 2, final time 
f = 20/32, At = h. 
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Figure 2: Average, minimum and maximum number of Newton iterations per- 
formed during the integration until final time. In (a) At was kept fixed and m 
varied, in (b) m = 2 was kept fixed and At varied. 



4.1 Convergence of the global method and of Newton's 
method 

First we check the convergence of the method. We perform test for m ranging 
from 2 to 5, observing no appreciable difference in the convergence properties 
of the algorithm. In all tests we choose At = h. 

Figured] plots the l-i errors between the numerical solution at time t — 20/32 
and the exact solution (|38]l and shows that the method is first order convergent, 
as expected for this choice of time stepping procedure and also due to the 
presence of the singularity in the first derivative of the exact solution. The 
dashed line is a reference slope for first order schemes. We observe that the 
convergence is not significantly affected by the parameter m. 

Figure [2] plots the number of Newton iterations employed by the algorithm 
during the integration from t — to t — 20/32. We plot the average (circles), 
minimum and maximum (solid lines) number of Newton iterations per timestep. 
Taking At — h (Figure [2^) , we observe that the number of Newton iterations 
slowly decreases when N increases and that, for any given N it increases only 
very moderately when m increases. In the case m = 2 we also tried to vary the 
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First timestep with M=2 (black solid), M=3 (rod dashed) 
10° | 1 1 1 
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Figure 3: History of the convergence of the Newton iterations during the first 
timestep. Black solid lines correspond to m = 2 and red dashed ones to m = 3. 
We show the results for N ranging from 32 to 4096, with At = h: the behaviour 
under grid refinement is indicated by the thin arrow. 

step size from At = h/10 to At = 5h. The results are reported in Figure [2b, 
showing that the number of Newton iterations grows when taking larger At in 

. The larger variability (for fixed N) and the irregular behaviour of the mean 
value when increasing TV in the case At = 5h preludes to the loss of convergence 
that we observe if At is taken even larger. 

Next we verify the convergence of the Newton's method. In Figure [3] we 
plot the Newton's error estimate ||u 1,fe+1 — u 1 -' c ||/||u 1 ' fe || obtained when com- 
puting the first timestep u 1 . We compare different number of grid points 
(N = 32, 64, . . . , 4096) as indicated by the thin arrow and two values for the 
exponent m appearing in (J3TJ) . 

We emphasise that as prescribed in Proposition 12.71 the choice of At = h is 
acceptable for the convergence both of the global numerical scheme and for the 
convergence of the Newton procedure. 

4.2 Solution of the linear system 

This section is devoted to computational proposals for the solution of a linear 
system where the coefficient matrix is the Jacobian in (|7J|, which is required at 
every step of the Newton procedure. For all the tests, we set m = 2, final time 
t = 20/32, At = h, and we let N be equal to 32, 64, ... , 1024 for checking the 
optimality of the proposed best solvers. 

As already stressed in Rcmark l2.ll the matrix is (weakly) non-symmetric so 
we start by considering the use of preconditioned GMRES (PGMRES). 

4.2.1 GMRES 

In Figure |4^, we plot the average (circles), minimum and maximum (vertical 
lines) number of GMRES iterations performed during the integration until final 
time, at different spatial resolutions. A least square fit (dashed line) shows that 
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GMRES iterations 
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(a) 
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Figure 4: Average, minimum and maximum number of GMRES iterations (a) 
and preconditioned GMRES iterations (b) performed during the integration 
until final time. The dashed line in panel (a) is the least square fit. 



the number of iterations grows as jV 5320 . This fact is in complete accordance 
with the analysis of Subsection 13.21 and in particular with equation (|34p . 

In Remark 13.51 we proved that Yat(u) is negligible with respect to the sym- 
metric positive definite term Ajv(u). Accordingly, in Remark 13.61 the use of 
Xn(u) as preconditioncr for F'(u) was analysed and it was shown to provide a 
strong spectral clustering of the preconditioned matrix at 1 and therefore we ex- 
pect a number of iterations not depending on the size N of the matrix as in ()35j) : 
this fact is observed in practice and indeed the iteration count of the PGMRES 
is almost constant, with average value equal to 6 iterations (see Figure 0})). 

At this point we are left with the problem of solving efficiently a generic 
linear system with coefficient matrix Xjv(u), which is a regularised version of a 
weighted Laplacian (i.e., by re-scaling, it is a shift of — Ld{u) by h 2 /At times 
the identity). A standard V-cycle is thus optimally convergent since Xn(u) is 
slightly better conditioned than a standard Laplacian. 

4.2.2 CG 

Since the non-symmetric part of F'(u) is negligible, we can try directly the 
solution of the whole system by using techniques such as the preconditioned 
CG (PCG) or the multigrid method which in theory should suffer from the loss 
of symmetry in the linear system. 

In Figure [5^, we plot the average (circles), minimum and maximum (vertical 
lines) number of CG iterations performed during the integration until final time, 
at different spatial resolutions. A least square fit shows that the number of 
iterations grows as jV ' 5 . As previously observed in connection with the 
GMRES method, the number of iterations is again essentially proportional to 
\/N, which agrees with the discussion in Subsection 13.21 

However, the number of iterations for fine grids is a lot higher that the ones 
with GMRES (up to 950 instead of 160 with a grid of 1024 points) so the latter 
has to be preferred. 
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Figure 5: Average, minimum and maximum number of CG iterations (a) and 
PCG iterations (b), performed during the integration until final time. The 
dashed line in panel (a) is the least square fit. 
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Figure 6: Average, minimum and maximum number of MGM iterations per- 
formed during the integration until final time. 



In a similar way, we consider A/v(u) as preconditioner in the PCG method. 
Results are shown in Figure EJd. The number of iteration is again essentially 
constant with respect to N, but also in the preconditioned version, the GMRES 
is slightly better since the number of PCG iterations, 8 or 9, is higher that the 
number of PGMRES iterations which was equal to 6. Furthermore we have a 
higher variance in the number of iterations, due to the weak non-symmetry of 
whole matrix. 

4.2.3 MGM 

We test the optimality of MGM, as discussed in Section [3~2l We apply a single 
recursive call, that is the classical T^-cycle procedure. As smoother, we use a sin- 
gle Jacobi step with damping factor equal to 2/3. We observe mesh independent 
behaviour with 10 or 11 iterations (see Figure |6]). 

We also tried other more sophisticated multi-iterative approaches by adding 
one step of post-smoother with Gauss-Seidel or standard Jacobi: the number 
of iterations drops to 6, but the cost per iteration is almost doubled, so that 
we do not observe a real advantage. The use of one step of CG or one step 
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Figure 7: Average, minimum and maximum number of PGMRES (a) and PCG 
(b) iterations performed during the integration until final time, when using one 
MGM V-cycle as preconditioner. 



of GMRES as post-smoother is not better while one step of PGMRES with 
preconditioner equal to Xn(u) reduces the number of iterations, but not enough 
compared with the cost of the solver needed for handling a generic system with 
the preconditioner as coefficient matrix. 

4.2.4 Krylov methods with MGM as preconditioner 

The previous experiments confirm that the MGM is an excellent solver for our 
linear system. Often this method is also applied as preconditioner in a Krylov 
method instead of employing it as a solver. In other words, as preconditioning 
step, we perform a single V-cycle iteration, with the same coefficient matrix and 
where the datum is the residual vector at the current iteration. 

With the use of such very cheap MGM preconditioning, the PGMRES con- 
verges within 4 or 5 iterations independently of the size of the involved matrices 
(see Figure Uh-)- Comparing with the GMRES method preconditioned with the 
symmetric part of F'(u) considered in l4.2.Il and Figure0b, the present precondi- 
tioning strategy is not only computationally cheaper, but it is also more effective 
since it achieves a stronger reduction of the number of GMRES iterations. 

Analogously, the application of the same MGM preconditioning in the PCG 
method leads to a convergence within 7 or 8 iterations, again independently of 
the system sizes (see Figure [7b). 

In conclusion, V-cycle preconditioning in connection with GMRES has to 
be preferred, taking into account the simplicity, the robustness (less variance 
in the iteration count), and the number of iterations. Indeed, due to the small 
iteration count, also the memory requirement does not pose any difficulty, since 
the number of vectors that have to be stored in the GMRES process is very 
reasonable. 
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5 2D generalization 



In this section we describe a straightforward 2D generalization of the numerical 
approach studied in the previous part of the paper. To this end, we consider 
a rectangular domain = [ao,ai] x [60,61] C R 2 and the grid points Xij = 
(ao + ih, bo + jk). For simplicity and without loss of generality, we also assume 
that the region f2 is square and choose identical discretization steps in the two 
directions (i.e. h — k), so that using N points per direction we have h = k = 
(a 1 -a )/(N + l) = (bi-b )/(N+l). The grid is thus composed of the (N + 2) 2 
points Xi j for i and j ranging from to iV+1. We denote with Uij the numerical 
value approximating u(xij). Of course, as in the one-dimensional case the use 
of Dirichlet boundary conditions allows to reduce to gridding to the N 2 internal 
points. 

In this setting, we generalize the finite difference discretization ([3]) of the 
differential operator as follows: 

9 ( t — . / ^du\ d / . .du 

Di+l/2,jUi+l,j — (D i+1 /2,j + Di_i/2,j) Uij + Di_i/2jUi-l t j 
Di,j+l/2Ui,j+l — (Aj + 1/2 + Aj-l/2) Ui,j + A,j-l/2«ij-l 



h 2 

where we denoted 



•o(l), (39) 



n _ Di+i,j + Di,j n _ Ih., ■ 1 + Dj,j 

U i+l/2,j - 2 ^,3+1/2 - ^ ' 

In order to write in matrix form the approximated differential operator above, 
we must choose an ordering of the unknowns Uij, arranging them into a vector 
u and approximate 

[V • (D(u)Vu)( Xl . 3 )}l J=1 ~ ^L D(n) u. 

The positions of the nonzero entries of the matrix L D ^ of course depend on 
the chosen ordering, so here we keep a double-index notation for the elements 
of u and of the matrix entries. Therefore, following (|39|) . L D ^ has entries 

[ifl(n)]jj = fii,l$j,m ( — A+1/2J — A-1/2J — Aj + 1/2 — Aj-l/2) 
+ Sl,i+lfim,jDi+l/2,j + fil,i-lfim,jDi_i/2J 

+ AiJ+l-DjJ+1/2 + fil,i3m,j-lDij_i/2 

on the (i,j) th row and (l,m) th column. The actual sparsity pattern of the 
resulting matrix thus depends on the ordering of the unknowns u^j] with the 
usual lexicographic ordering that has u^j in the (i + N(j — 1)) position of 
the vector u, one may have, as in the case of the standard Laplacian operator, 
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nonzero entries only on the main diagonal, on the 1 st and N th upper and lower 
diagonals. 

Each timestep with the Crandall Liggett formula ^ thus requires the solu- 
tion of the nonlinear equation defined by 

F(u) = u - ^r L D(u) u - u™" 1 . 

As in the one-dimensional case we propose to approximate the solution of the 
nonlinear equation with the Newton's method; an analysis similar to that of 
Theorem 12 . 71 can be carried out in the new 2D context. The Jacobian of F(u) 
is 

F'(u)=I-—L D(u) -—Y(u) (40) 



where 



h 2 h 2 

i l.m 



Y i,j ( u ) = L dUrs Ul > m ( 41 ) 



l.m 

with the double-index notation as above. 

A tedious but straightforward computation yields 

+ 2 D 'i+l,]^r,i+l5 s , 3 {U i+ Xj - Uij) + -D'^jS^i-lSs^ (Ui-l,j - U^j) 

+ -D' i j+1 S r ^5 s , J+ i {uij+i - mj) + -Dl j _ 1 S„S s j-i (iM,j-i - Ui,j) 

for the generic entry of Y(u). (The obvious changes must be taken into account 
to implement the boundary conditions, e.g. either eliminating the unknowns for 
the points on the Dirichlet boundary or the unknowns on suitably chosen ghost 
points outside the Neumann boundary.) 

As in the one-dimensional case, the matrix Y(u) may be written as ^(u) = 
T(u)D'(u) where D'(u) is a diagonal matrix with entries equal to D'(uij) and 
in smooth regions of the solution, the nonzero entries of Y(u) are 0(h 2 ) on 
the main diagonal and 0(h) outside. Moreover, entries of Y(u) expected to 
be nonzero may in fact be null because the approximate solution is locally flat 
in a neighbourhood or because some of the D'(uk) may be null. When using 
the natural ordering of the unknowns described above, the sparsity patterns of 
Y(u) and F'(u) for the Barenblatt solution are illustrated in figured] The gaps 
along the diagonals of Y(u) correspond to the regions where the approximate 
solution u is flat. 

We perform ed our tests with the two-dimensional Barenblatt solution (see 
Vazquez , 20071 ) with exponent m = 4 on grids of size N x N for N ranging from 



32 to 1024. First of all we note that the number of Newton iterations required at 
each timestep is almost independent of N and is (on average) 4 when At = 0.5h, 
4.5 when At = h and 6.5 when At = 2h (see Figure 
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Sparsity pattern of Y 



Sparsity pattern of J 




V 



(b) 



\ 

\ 



Figure 8: Sparsity pattern of Y(u) (a) and F'(u) (b) on a 10 x 10 grid with the 
unknowns in lexicographic ordering. 




Figure 9: Average, minimum and maximum number of Newton iterations per- 
formed during the integration until final time. The 3 data series for each ./V 
have been slightly shifted for clarity. 




Figure 10: Number of GMRES iterations at different grid sizes, in 2D. (a) 
without preconditioning, (b) with V-cycle preconditioner. On the right, the 3 
data series for each N have been slightly shifted for clarity. 
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We point out that the results of Section |3~T1 generalize to the two-dimensional 
case and thus we perform numerical tests using a multigrid iteration as pre- 
conditioner for PGMRES which in Section @] provided best results in the one 
dimensional case. 

In particular we employ a single V-cycle iteration, with a Galerkin approach 
using the bilinear interpolation as prolongation operator and one step of red- 
black Gauss-Seidel as pre-smoother. In Figure [10] we plot the mean (symbols) 
and minimum- maximum (solid lines) number of GMRES iterations needed at 
different spatial resolutions. Different colours correspond to different choices 
of At, namely At = h/2 (blue crosses), At — h (black circles) and At = 
2h (red diamonds). The left panel shows that, without preconditioning, the 
number of GMRES iterations grows with the grid size: least square fits yield 
the approximations TV - 5165 , TV - 5435 and jV 5702 respectively for the number of 
GMRES iterations on an N x N grid with the three choices of At mentioned 
above. For homogeneity, the results for N — 1024 are not reported in the graph, 
since they require the restarted GMRES method or a parallel implementation, 
due to memory limitations when run on a PC with 8Mb of RAM. 

Figure [TOb clearly demonstrates the optimality of the preconditioning strat- 
egy adopted, with the number of iterations being in the narrow range 5-10 when 
TV ranges from 32 to 1024 and with all the three choices of the time step and 
with the average number of iterations being always between 5 and 7. We note 
in passing that we also employed damped Jacobi as a smoother with analogous 
results on the optimality, but observing a slightly higher number of iterations 
(8-11 on average). 



6 Conclusions and future developments 

The novel contribution of this paper relies in the proposal of a fully implicit 
numerical method for dealing with nonlinear degenerate parabolic equations, in 
its convergence and stability analysis, and in the study of the related computa- 
tional cost. Indeed the nonlinear nature of the underlying mathematical model 
requires the application of a fixed point scheme. We identified the classical New- 
ton method in which, at every step, the solution of a large, locally structured, 
linear system has been handled by using specialised iterative or multi-iterative 
solvers. In particular, we provide a spectral analysis of the relevant matri- 
ces which has been crucial for identifying appropriate preconditioned Krylov 
methods with efficient V-cycle preconditioners. Numerical experiments for the 
validation of our multi-facet analysis complement this contribution. 

Among the vast range of possible applications of degenerate parabolic equa- 
tions, we poin t out a recent one in the field of monument conservation in 
ISemplice et al. (2009), where an approximation technique derived from the one 



analysed here has been successfully employed in the forecast of marble deteri- 
oration on monuments. Having in mind the application to more complicated 
monument geometry, we will pursue the extension of the results of this paper 
to the case of finite element methods for the space discretization. 
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